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CONTRACTOR:  Massachusetts  Institute  of  Technology 


PRINCIPAL  INVESTIGATOR:  James  C.  Weaver 


CONTRACT  TITLE:  Electroporation:  Theory  of  Basic  Mechanisms 
INTRODUCTION 


The  objective  of  our  investigation  is  the  development  of  a  theory  of  the  mechanism  of  electroporation. 
Overall  we  seek  to  find  models  which  allow  theoretical  description  of  measurable  quantities.  This 
includes  the  more  specific  objective  of  first  providing  a  quantitative  description  of  key  features  of 
electrical  behavior,  and  subsequently  of  molecular  transport. 

I  1  ‘ 

'‘in  summary  form,  electroporation  (i)  is  now  believed  to  be  a  universal  cell  membrane  phenomenon, 
involving  both  the  lipid  bilayer  and  membrane  macromolecules,  and  is  therefore  fundamental  to  mem¬ 
brane  understanding,  (ii)  provides  a  general  method  for  introducing  molecules  into  cells,  or  releasing 
molecules  from  cells,  with  potentially  major  applications  in  science  and  technology,  and  yet  (iii)  its 
mechanism  is  poorly  understood.1'3  For  example,  no  previous  theory  actually  describes  electrical 
behavior  during  electroporation,  membrane  recovery,  or  the  amount  of  molectifar  transport. 


APPROACH  TO  A  THEORETICAL  MODEL 
With  this  in  mind,  our  specific  goals  have  been: 


V, 

t  iV.ftU.v,  w \  y  r;  '  * 
^  ^  v  U  V  OC  U  tv  v  .  r  r 


(1)  Extension  of  an  initial  theory  of  reversible  electrical  breakdown  to  one  with  more  solid  founda¬ 
tions,  i.c.  elimination  of  the  approximate  "switch  on"  criteria  of  pores,  and  elimination  of  the 
assumption  of  the  membrane  containing  so  many  pores  that  the  membrane  was  "saturated”  with 
pores.4  A  significant  advance  towards  this  goal  has  been  made,  and  is  presented  in  a  recently 
submitted  paper  (copy  appended). 


(2)  Achievement  of  a  theory  which  quantitatively  describes  the  transmembrane  potential,  U(t),  during 
irreversible  rupture,  such  that  a  unified  theory  of  both  REB  and  rupture  is  provided  by  one 
model.  Such  a  theory  should  yield  predictions  of  U(t)  which  can  be  compared  directly  with 
experiments,  a  basic  requirement  of  which  has  not  yet  been  achieved  by  other  theories  of  electro¬ 
poration.  This  has  recently  been  accomplished,  and  is  presented  in  the  appended  manuscript 


(3) 


Related  achievement  of  the  ability  to  quantitatively  describe  incomplete  REB,  i.e.  a  discharge  that 
stops  before  the  transmembrane  potential  reaches  zero.  This  has  also  been  achieved  (see 
appended  manuscript). 
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(4)  Extension  of  our  first,  successful  theory  of  the  reversible  electrical  breakdown  of  electroporation 
to  include  metastable  pores  associated  with  a  "foot-in-the-door"  mechanism,  i.e.  a  moleculey^ore 
interaction  wherein  the  presence  of  a  macromolecule  partially  inserted  within  a  pore  prevents  the 
pore  from  shrinking  and  then  vanishing.  This  approach,  based  on  the  approximation  of  one¬ 
dimensional  diffusion  through  the  long  lifetime  pore,  is  a  candidate  mechanism  for  molecular 
transport  associated  with  electroporation. 

(5)  Further  extension  to  include  a  pore-membrane  macromolecule  interaction  in  the  theory.  Mem¬ 
brane  channel  proteins  are  prime  candidates  for  such  interactions,  and  may  provide  nucleation 
sites  for  the  long  lifetime  pores  which  are  believed  to  occur.  Such  mctastablc  pores  may  have 
significantly  longer  lifetimes  than  "ordinary"  transient  aqueous  pores,  and  may  exhibit  a  strong 
temperature  dependence  of  the  lifetime.  This  candidate  mechanism  may  also  be  relevant  to 
describing  molecular  transport. 

(6)  As  partially  noted  above,  carry  out  further  work  towards  the  development  of  an  extended  theory 
which  predicts  both  electrical  behavior  and  the  amount  of  transmembrane  transport  of  molecules. 
Here  the  involvement  of  membrane  proteins,  particularly  channel  forming  proteins,  is  believed  to 
be  important.  Such  a  theory  should  include  hindered  diffusion  as  a  primary  mechanism,  particu¬ 
larly  for  long  times  involving  persistent  metastable  pores.  However  electrokindic  mechanisms 
such  as  electrophoresis  and  clcctroosmosis  which  may  also  operate  during  the  time  the  transmem¬ 
brane  potential  is  non-zero.  The  general  goal  is  to  predict  the  number  of  molecules  which  move 
across  a  cell  membrane,  and  also  (because  of  the  fundamentally  statistical  nature  of  the  theory, 
and  the  statistical  orientation  of  non-spherical  alls)  the  distribution  of  transport  within  a  cell 
population.5 

Throughout  we  have  sought  quantitative  estimates  which  can  be  compared  to  the  results  of  experiments 

by  ourselves  and  others. 


SUMMARY 

We  have  completed  our  initial  extension  of  a  transient  aqueous  pore  theory  of  electroporation.  This 
improved  model  yields  descriptions  of  four  key  aspects  of  complex  electrical  behavior 

•  Reversible  electrical  breakdown  (REB)  leading  to  complete  membrane  discharge 

•  Incomplete  REB  (discharge  halts  at  U  >  0) 

•  Rupture  (mechanical)  with  its  characteristic  slow,  sigmoidal  electrical  discharge 

•  Membrane  charging  without  dramatic  behavior  at  small  U 

The  extended  theory  quantitatively  describes  this  complex  set  of  behavior.  This  improved  version 
eliminates  the  use  of  sharp  conduction  criteria  for  pores,  and  instead  uses  a  more  realistic  continuous 
conduction  in  which  an  estimate  of  the  Bom  energy  is  used  with  a  Boltzmann  factor  to  describe  the 
reduced  conductivity  of  a  pore  in  a  low  (compared  to  water)  dielectric  constant  membrane.  This  ver¬ 
sion  also  eliminates  the  use  of  an  assumption  concerning  the  number  of  pores  present  in  a  membrane 
at  equilibrium,  and  instead  utilizes  creation  and  destruction  rates  for  pores.  This  version  further 
assumes  a  realistic  minimum  pore  size,  based  on  molecular  sizes.  Finally,  this  version  provides  the 
first  unified  and  quantitative  description  of  the  several  dramatic  electroporation-related  phenomena 
listed  above. 
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We  are  now  continuing  to  make  progress  towards  combining  the  improved  version  (which  describes 
the  very  early  electrical  events)  with  a  model  for  pore-membrane  protein  interactions  (expected  to  be 
relevant  to  the  long-lived  high  permeability  state).  Although  still  fraught  with  difficulties,  we  antici¬ 
pate  that  this  approach  will  lead  to  a  reasonable  description  of  cell  membrane  electroporation,  wherein 
first  dramatic  electrical  events  occur,  followed  by  two  phases  of  membrane  recovery  (fast,  complete 
resealing  of  most  pores,  and  slow,  metastable  pore-protein  complexes  which  recover,  probably  through 
a  thermally  activated  process).  Finally,  we  expect  that  mass  transport  through  both  the  fast  recovering 
pores  and  slowly  recovering  pores  will  be  significant,  with  a  combination  of  hindered  diffusion  and 
electrokinctic  transport  being  significant  for  the  fast  recoverying  pores,  and  primarily  hindered  diffusion 
being  important  through  the  slowly  recovering  pore-protein  complexes.  Throughout,  we  will  continue 
to  emphasize  theoretical  development  which  can  lead  to  the  prediction  of  experimental  results,  thereby 
allowing  direct  comparisons. 


COPY  OF  SUBMITTED  MANUSCRIPT 

Because  of  the  complexity  and  length  (47  pages)  of  the  calculations,  text  and  computed  graphs  which 
show  how  one  can  quantitatively  describe  the  electrical  behavior  associated  with  electroporation  in  an 
artificial  planar  bilayer  membrane,  we  have  not  attempted  to  summarize  the  calculation  and  its  results, 
as  this  could  be  misleading.  Instead  we  have  provided  a  complete  copy  of  the  manuscript: 

Barnett,  A.  and  J.  C.  Weaver  "Electroporation:  A  Unified,  Quantitative  Theory  of  Reversible 
Electrical  Breakdown  and  Rupture"  (submitted). 

Figures  5  through  12  provide  illustrations  of  the  behavior  of  the  model,  and  demonstrate  that  much  of 
the  complex  behavior  observed  in  experiments  is  actually  described  by  the  model. 


References 

1.  Zimmermann,  U.,  ‘‘Electrical  Breakdown,  Electropermeabilization  and  Electrofusion,”  Rev.  Phy¬ 
siol.  Biochem.  Pharmacol.,  vol.  105,  pp.  175  -  256,  1986. 

2.  Potter,  H.,  ‘‘Electroporation  in  Biology:  Methods,  Applications,  and  Instrumentation,”  Analyt. 
Biochem.,  vol.  174,  pp.  361  -  373,  1988. 

3.  Neumann,  E.,  A.  Sowers,  and  C.  (Eds)  Jordan,  Electroporation  and  Electrofusion  in  Cell  Biology, 
Plenum,  New  York,  1989. 


4.  Weaver,  J.  C.  and  K.  T.  Powell,  ‘‘Theory  of  Electroporation,”  in  Electroporation  and  Electrofu¬ 
sion  in  Cell  Biology,  ed.  E.  Neumann,  A.  Sowers  and  C.  Jordan,  pp.  Ill  -  126,  Plenum,  1989. 


'ffc 

Electroporation:  A  Unified,  Quantitative 
Theory  of  Reversible  Electrical  Breakdown  and  Mechanical  Rupture 
in  Artificial  Planar  Bilayer  Membranes 

Alan  Barnett  and  James  C.  Weaver 
Harvard-MIT  Division  of  Health  Sciences  and  Technology 
Cambridge,  MA  02139,  U.S.A. 


Aoril  5.  1989 


Page  2 


ABSTRACT 


We  present  a  quantitative  theory  of  electroporation  of  artificial  lipid  bilayer  membranes.  Assum¬ 
ing  that  aqueous  pores  cause  electroporation,  we  describe  the  pore  population  of  the  membrane  by  the 
density  function  n(r,r),  where  n{r,t)dt  is  the  number  of  pores  with  radius  between  r  and  r  +  dr  at 
time  r.  To  write  a  set  of  differential  equations  for  the  evolution  of  n(r,t),  we  assume  that  there  is  a 
minimum  pore  size  rmin,  that  pores  of  radius  rmin  are  created  and  destroyed  by  thermal  fluctuations, 

U 2 

and  that  the  pore  creation  rate  is  proportional  to  exp(-  where  U  is  the  membrane  voltage,  we 


derive  a  relation  between  pore  size  and  pore  conductance,  we  use  the  expression  for  the  pore  energy 
previously  derived  by  Pastushenko  and  Chizmadzhev,  and  we  include  a  model  of  the  external  circuit, 
equation  numerically.  We  solve  the  equations  numerically  and  compare  the  solutions  to  the  results  of 
charge  pulse  experiments. 

In  a  charge-pulse  experiment  a  membrane  suffers  one  of  four  possible  fates:  (1)  a  slight  increase 
in  electrical  conductance,  (2)  mechanical  rupture,  (3)  partial  reversible  electrical  breakdown,  resulting 
in  incomplete  discharge  of  the  membrane,  or  (4)  reversible  electrical  breakdown  (REB),  resulting  in 
complete  discharge  of  the  membrane.  In  agreement  with  experiment,  our  theory  describes  these  four 
fates  and  predicts  that  the  fate  in  any  particular  experiment  is  determined  by  the  properties  of  the  mem¬ 
brane  and  the  duration  and  amplitude  of  the  charging  pulse. 


INTRODUCTION  AND  BACKGROUND 

Electroporation  is  a  set  of  related  phenomena,  caused  by  the  formation  of  aqueous  pores,  that  are 
observed  in  both  natural  and  artificial  bilayer  membranes  in  response  to  a  large  applied  electric  field. 

In  planar  artificial  membranes  these  phenomena  include  reversible  electrical  breakdown  (REB)  and 
membrane  rupture1,  while  in  cells  REB  is  followed  by  a  transient  high  permeability  state.2'5  Much  of 
the  recent  interest  in  electroporation  relates  not  to  electrical  behavior,  but  rather  to  the  significant  fluxes 
of  molecules  across  cell  membranes  that  occur  during  the  high  permeability  state.6  Most  applications  of 
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electroporation  have  involved  the  introduction  of  genetic  material  into  cells,  but  the  large  transient 
molecular  fluxes  have  much  wider  applications.  Although  such  molecular  transport  itself  has  some¬ 
times  been  termed  "electroporation”7,  the  transient  molecular  flux  is  only  one  consequence  of  electri¬ 
cally  generated  large  pores;  wc  therefore  use  the  term  "electroporation"  to  refer  to  all  pore  phenomena 
associated  with  a  large  transmembrane  potential  difference. 

Although  the  topology  of  cells  differs  from  that  of  planar  membranes,  many  features  of  REB  in 
planar  membranes  and  in  cell  membranes  are  the  same.8  The  planar  geometry,  which  permits  easy 
access  to  both  sides  of  the  membrane,  and  the  simplicity  of  bilayer  membranes  that  lack  proteins  com¬ 
bine  to  make  artificial  membranes  a  good  model  system.  The  insight  gained  from  the  study  of  artificial 
membranes  will  be  a  good  foundation  on  which  to  build  a  theory  of  electroporation  in  biological  mem¬ 
branes.  The  remainder  of  this  paper  deals  exclusively  with  artificial  bilavcr  membranes. 

First  we  discuss  two  types  of  experiments  used  for  studying  electroporation  of  artificial  planar 
bilayer  membranes,  then  we  describe  phenomena  observed  in  the  experiments,  and  finally  we  derive 
and  discuss  a  theory  that  can  explain  many  observations  by  providing  a  complete  description  of  the 
observable  electrical  behavior  of  artificial  bilayer  membranes,  including  the  trans-membrane  potential 
U(t)  and  the  membrane  conductance  G(t). 

ELECTROPORATION  EXPERIMENTS:  VOLTAGE  CLAMP  AND  CHARGE  PULSE 

There  are  two  major  types  of  electroporation  experiments:  voltage  clamp  and  charge  pulse.  Both 
types  of  experiment  use  apparatus  similar  to  that  shown  in  Figure  I.  The  apparatus  consists  of  a 
vessel  filled  with  an  electrolytic  solution.1- E'1-  A  planar  lipid  bilayer  membrane  spans  the  vessel,  divid¬ 
ing  it  into  two  compartments.  In  each  compartment  is  a  planar  electrode  oriented  parallel  to  the  mem¬ 
brane.  When  the  electrodes  are  at  different  potentials,  the  electric  field  is  perpendicular  to  the  mem¬ 
brane  .  An  external  circuit  applies  a  signal  to  the  electrodes  and  measures  the  response  of  the  system. 
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An  equivalent  circuit  for  a  voltage  clamp  experiment  is  shown  in  Figure  2  and  for  a  charge  pulse 
experiment  is  shown  in  Figure  3.  We  model  the  membrane  as  a  capacitor  C  in  parallel  with  a  resistor 

— - — .  The  membrane  is  in  series  with  the  resistor  Re  which  represents  the  electrolyte  and  the  elec- 

<7(0 

nodes  (we  neglect  the  capacitance  of  the  electrodes,  which  is  small).  The  external  circuit  consists  of  a 
voltage  source  with  internal  resistance  RT  in  series  with  an  ammeter  for  a  voltage  clamp  experiment 
and  of  a  current  source  with  internal  resistance  RN  in  parallel  with  a  voltmeter  for  a  charge  pulse 
experiment. 

In  a  voltage  clamp  experiment,  a  voltage  waveform  (7,  (0.  typically  a  scries  of  step  functions,  is 
applied  to  the  electrodes  and  the  current  Im ,  the  sum  of  the  charging  current  to  the  "membrane  capaci¬ 
tor"  and  the  leakage  current  through  the  "membrane  resistor”is  measured.  The  quantity  of  interest  is 
the  membrane  conductance  G  (r).  When  the  applied  voltage  is  not  changing  G  (r )  is  related  to  Ut  (r ) 
and  to  the  measured  current  lm(t)  by  the  equation 

U t  (i )  —  Im  (Re+Rt) 

G  (0  = - — - — — T-~—  (1) 

!m  +  (7?£  +  7?r)C  d~ 

Typically,  the  amplitude  of  the  applied  voltage  is  of  the  order  of  0.2  v  -  0.7  v,  and  the  duration  of  the 
pulse  is  of  the  order  of  milliseconds  or  longer.8 


In  a  charge  pulse  experiment  a  current  waveform  /,(r),  typically  a  square  wave,  is  applied  to  the 
electrodes  and  at  the  end  of  the  pulse  the  switch  is  opened.1  After  the  pulse,  the  membrane  conduc¬ 
tance  G  (: )  is  related  to  the  measured  voltage  Um  by  the  equation 


G(t)  =  -C 


d(\nUn) 

dt 


(2) 


Since  a  membrane  undergoing  REB  can  discharge  in  less  than  a  microsecond,  the  time  scale  of 
charge  pulse  experiments  is  much  shoncr  than  the  time  scale  of  voltage  clamp  experiments.  Although 
the  theory  we  develop  is  applicable  to  the  behavior  of  membranes  on  all  time  scales  of  interest,  in  this 
paper  we  confine  ourselves  to  the  discussion  of  charge  pulse  experiments  and  shon  time  scales. 
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MEMBRANE  FATES  IN  CHARGE  PULSE  EXPERIMENTS 

Benz  et  a/1  have  observed  that  a  membrane  in  a  charge  pulse  (pulse  length  =  0.4  ps)  experiment 
suffers  one  of  four  possible  fates: 

(1)  Slight  increase  in  conductivity 

(2)  Irreversible  rupture 

(3)  Incomplete  reversible  electrical  breakdown,  or 

(4)  Reversible  electrical  breakdown. 

In  case  4,  the  membrane  conductance  increases  by  up  to  eight  orders  of  magnitude,  causing  the  mem¬ 
brane  to  discharge  completely  in  less  than  1  (is.  The  membrane  then  recovers  slowly  to  its  original 
state.  In  case  3,  the  membrane  conductance  increases  by  several  orders  of  magnitude;  the  membrane 
then  begins  to  discharge  with  a  characteristic  time  of  about  1  ps,  but  the  membrane  recovers  before  the 
membrane  discharges  completely.  In  case  2,  the  membrane  conductance  increases  sharply  after  a  delay 
of  many  microseconds  and  never  recovers  its  original  properties  following  discharge.  In  case  1 .  the 
membrane  conductance  increases  slightly  during  the  pulse  and  slowly  returns  to  its  original  value.  The 
fate  of  the  membrane  in  any  given  experiment  is  determined  by  the  properties  of  the  membrane  and  the 
duration  and  magnitude  of  the  applied  pulse.  A  successful  electroporation  theory  must  predict  the  out¬ 
come  of  any  given  experiment  and  explain  the  differences  between  the  four  fates  in  terms  of  a  reason¬ 
able  physical  model. 

In  the  present  theory  wc  hypothesize  that  a  membrane  contains  a  population  of  aqueous  conduct¬ 
ing  pores  that  arc  created  and  destroyed  by  thermal  fluctuations.  Our  "standard  membrane"  in  equili¬ 
brium  at  room  temperature  and  with  membrane  voltage  U  =  0  contains  about  seven  pores.  We  will 
show  that  our  theory  gives  the  following  explanation  for  the  four  fates: 

If  a  small  amplitude  pulse  is  applied  to  the  membrane  experiment,  the  membrane  charges  with  a 
characteristic  RC  time  constant.  The  membrane  resistance  does  not  increase  significantly,  as  the  pore 
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creation  rate  will  not  change  much  from  its  value  for  U  =  0.  If  the  pulse  amplitude  and  duration  are 
large  enough  to  cause  rupture,  the  electric  field  across  the  membrane  increases  the  pore  creation  rate 
and  causes  the  resulting  pores  to  increase  in  size.  The  membrane  conductance  increases,  but  not 
enough  to  cause  the  membrane  to  discharge  rapidly.  Eventually,  one  or  more  pores  become  large 
enough  to  become  unstable.  As  the  large  pore  grows,  the  membrane  conductance  increases  and  the 
membrane  discharges,  but  the  unstable  pore  continues  to  grow  until  the  membrane  ruptures.  If  the 
pulse  amplitude  and  duration  arc  large  enough  to  cause  REB,  the  enhancement  of  the  membrane  con¬ 
ductance  is  sufficient  to  discharge  the  membrane  before  any  pores  become  unstable;  the  pores  then 
shrink  in  size  and  their  number  decreases,  and  the  membrane  returns  to  its  initial  state  undamaged. 
For  partial  REB,  the  discharge  is  incomplete;  the  membrane  conductance  increases  sufficiently  rapidly 
that  the  membrane  remains  charged  to  .  1  v  or  so. 


MATHEMATICAL  FORMULATION 

We  now  present  a  derivation  of  our  theory,  which  is  an  extension  and  improvement  of  previous 
work.1-'15  Consider  a  membrane  in  an  electroporation  experiment,  such  as  shown  in  Figure  1.  Let 
n{r ,  i)dr  be  number  of  pores  in  the  membrane  with  radius  between  r  and  r  +  dr  at  time  r.  The  den¬ 
sity  function  n(r,  r)  obeys  Smoluchowski 's  equation16 


on 

dr 


'd2n  ,  d 

n  d£E 

dr:  "  dr 

» 

kT  dr 

l  J  J 

(3) 


where  Dp ,  the  diffusion  constant  for  the  pore  radius,  is  independent  of  r 17.  k  is  Boltzmann’s  constant, 
T  is  the  absolute  temperature,  and  A E  is  the  pore  "energy",  a  function  of  r  and  U  with  dimensions  of 

3A£ 

energy  that  has  the  property  that  -(— - — )v  is  the  effective  force  acting  to  increase  the  pore  size. 

or 

Equation  (3)  is  valid  for  rmin  S  r  Srm]lx.  We  discuss  the  question  of  the  appropriate  values  of  rmm 
and  rmt,  in  the  section  on  boundary  conditions. 
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Equation  (3)  can  be  derived  from  the  assumption  that  j,  the  flux  density  of  pores  in  radius  space, 
is  given  by 


j 


bn  _n_ 
dr  kT 


bAE 

dr 


(4) 


The  first  term  in  the  right  member  of  Equation  (4)  is  the  "diffusion  flux  density"  due  to  random 
changes  in  the  size  of  the  pores  due  to  thermal  fluctuations.  The  second  term  is  the  "drift  flux  density" 
due  to  the  action  of  the  radial  force  acting  to  change  the  pore  size.  If  pores  are  neither  created  nor  des¬ 
troyed  (except  at  the  "boundaries”  rmin  and  rm„),  j  obeys  the  continuity  equation 


—  +  -^  =  0  (5) 

dt  dr 

Use  of  Equation  (4)  to  eliminate  j  in  Equation  (3)  yields  Equation  (1). 

To  use  Equation  (3)  to  describe  a  membrane  we  must  specify  the  functional  dependence  of  A E 
on  U  and  r,  apply  appropriate  initial  and  boundary'  conditions,  and  write  an  additional  equation  which 
relates  U(t)  to  n(r.t).  First  we  discuss  the  function  A E. 


THE  PORE  "ENERGY"  AE 

The  dynamics  of  the  pore  population  is  controlled  by  the  physical  forces  that  act  to  change  pore 
radius.  These  forces  are  of  two  types:  (I)  forces  due  to  thermal  fluctuations  that  change  rapidly  and 
randomly  in  time,  and  (II)  forces  that  are  functions  of  the  local  electric  field  and  the  mechanical 
configuration  of  the  pore.  A  complete  formulation  of  the  problem  would  require  knowledge  of  the 
stress  in  the  membrane  for  any  given  configuration.  For  a  given  membrane  and  pore  shape,  one  would 
solve  Maxwell  s  equations  for  the  electric  field,  and  use  the  Maxwell  stress  tensor  to  compute  the 
forces  of  electrical  origin.  Tnc  equilibrium  shape  of  the  membrane  and  pore  could  then  be  calculated 
by  balancing  the  mechanical  and  electrical  stresses.  One  would  then  have  to  repeat  the  calculation. 
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taking  into  account  thermal  fluctuations  of  the  pressure  in  the  bathing  solution,  the  stress  in  the  mem¬ 
brane,  and  the  electric  field.  This  approach  is  prohibitively  difficult.  To  make  the  problem  tractable, 
wc  follow  Abidor  ct  a/18  and  assume  that  each  pore  is  cylindrical,  and  is  therefore  completely 
described  by  one  parameter,  its  radius. 

We  write  A E  as  the  sum  of  the  mechanical  contribution  A £w,  the  electrical  contribution  A Ee, 
and  an  arbitrary  constant  A£0. 


A  E  —  AEm  +  A££  +  AE  Q  (6) 

A  simple  and  physically  reasonable  form  for  A EM  is19,20 

A E^  =  2rryr  -  rrl>2  (7) 

where  T  is  the  surface  tension  of  the  membrane-water  interface,  and  y  is  the  edge  energy  of  the  pore. 

To  compute  the  electrical  contribution  to  the  force,  wc  must  solve  the  electrostatics  problem. 

This  is  a  difficult  task,  and  no  adequate  formulation  of  the  problem  has  yet  appeared  in  the  literature. 
The  primary  difficulty  involves  computing  the  relation  between  the  current  density  and  the  electric  field 
in  the  aqueous  phase  near  the  membrane  and  inside  the  pores.  We  follow  Pastushcnko  and  Chizmad- 
zhev21  and  assume  that  the  system  can  be  described  by  Ohm’s  law,  and  that  the  electrical  conductivity 
c  is  constant  in  each  of  three  different  regions; 


c 


°<  in  the  bulk  electrolyte 
•j  c p  inside  a  pore 

0  inside  the  lipid 


(8) 


The  conductivity  cf  is  related  to  the  concentrations  C,  and  mobilities  t\i  of  the  ions  in  solution  by  the 


formula 


i 


(9) 
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where  t  -  1.6  x  10-19  coulomb,  z,  is  the  charge  of  the  i  -th  type  of  ion,  and  the  sum  is  over  all  of  the 
different  types  of  ions  in  the  solution.  The  conductivity  inside  the  pore  is  given  by 


op  =  £(z,e)2n,Cl//4exp 

4 

where  //,  ,  the  stcric  hindrance  factor,  is  a  function  of  the  pore  radius  and  of  the  radius  r,  of  ion  of 
type  i,  and  p,®  is  the  standard  chemical  potential  of  an  ion  of  type  i  inside  the  pore.  We  use  the  poly¬ 
nomial  approximation  for  //,  given  by  Rcnkin.22 


H(r) 


’ 

ft  » 

2' 

' 

»  « 

3 

*  « 

S' 

r, 

r4 

r, 

r,- 

1  - 

1-2.1 

+  2.09 

-0.95 

ft 

r 

» 

• 

r 

ft  4 

r 

r 

i  j 

ft 

(11) 


We  approximate  the  standard  chemical  potential  an  ion  in  the  pore  by  the  Bom  energy  of  a  point 
charge  on  the  axis  of  an  infinite  cylinder21 


=  ££  i- 

4  nz,r 


iL 


(12) 


where  e/  and  e»  arc  the  dielectric  constants  of  the  lipid  and  the  water,  respectively,  and  the  function  P 
has  a  maximum  value  of  0.25.23 

To  find  an  approximate  solution  to  the  electrostatics  problem,  we  start  with  the  assumptions  that 
the  pores  are  sufficiently  far  apart  that  their  mutual  interaction  is  negligible,  and  that  the  electric  field 
Ep  inside  a  pore  is  uniform  and  perpendicular  to  the  plane  of  the  membrane.  The  electrical  force  act¬ 
ing  to  enlarge  the  pore  can  be  expressed  in  terms  of  Ep.  The  electric  pressure  difference  at  the 
pore  edge  is 


kP'ltc 


(13) 


where  the  first  term  is  the  electrical  pressure  in  the  pore  and  the  second  is  the  electrical  pressure  in  the 
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adjacent  lipid.  The  electrical  contribution  to  the  pore  "energy"  A Ee  is  -1  times  the  work  done  to 
expand  the  pore  from  rmin  to  r ; 


A Ee  =  -  jAptltc  ( Inr'hdr ')  =  -n h(ew^tl)  j  E2r' dr'  (14) 

^*iin 

where  we  have  used  Equation  (13).  To  use  Equation  (14),  one  needs  to  find  an  expression  for  Ep  in 
terms  of  U .  Since  the  membrane  surface  is  not  an  equipotcntial  when  conducting  pores  arc  present, 
we  define  U  to  be  the  potential  difference  across  the  membrane  far  from  a  pore.  The  voltage  U  can 
be  written  as  the  sum  of  two  terms 


U  =  Up  +  Ut  (15) 

where  Up  is  the  voltage  drop  across  the  pore  itself  and  Us  is  the  voltage  drop  in  the  bulk  electrolyte 

near  the  ends  of  the  pore.  The  electric  field  Ep  is  related  to  Up  by 


E 

p  h 


(16) 


where  h  is  the  thickness  of  the  membrane.  The  voltage  drops  Up  and  Us  arc  related  to  /„ ,  the  current 


through  the  pore,  by 


and 


Up  =  Rplp 


U,  =  R,lp 


where  the  pore  resistance  Rp  and  the  spreading  resistance  Rs  24  are 


= 


Rs 


~r~cp 

1 

2oer 


Because  Rs  and  Rp  act  as  a  voltage  divider,  the  equation  relating  Ep  to  U  is 


Zp 


u  Rp 

h  Rs  +  Rp 


(17) 


(18) 

(19) 

(20) 


(21) 
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Using  Equation  (21),  we  can  now  rewrite  Equation  (14)  as 


A Ee  =  - 


*(6*  -t,)W 


J  a 2  r  dr 


(22) 


where 


a (r)  = 


1  + 


rcrqp(r) 

2hae 


-l 


(23) 


a  result  first  obtained  by  Pastushcnko  and  Chizmadzhev  21  Note  that  for  small  values  of  r,  a  =  1, 
while  for  large  r  both  a  and  the  electrical  force  tend  to  0  as  — . 


It  is  convenient  to  define  the  constant  A £0  such  that 


A £(rmin)  e  0  (24) 

The  function  A E(r,U)  is  completely  defined  by  Equations  (6)  *  (12)  and  (19)  -  (24).  Figure  4  is  a 
plot  of  A E  versus  r  for  various  values  of  U .  The  values  of  the  parameters  y.  T,  E/,  £»,,  h,  and  r, 
used  in  the  computation  for  Figure  4  are  given  in  Table  1. 

BOUNDARY  AND  INITIAL  CONDITIONS 

Since  Equation  (3)  is  a  parabolic  partial  differential  equation,  appropriate  boundary'  conditions 
consist  of  two  equations  of  the  form 


r  > 


J 


(25) 


one  valid  at  r  =  rmin  and  the  other  valid  at  r  =  rmix  for  all  times,  and  an  initial  condition  specify¬ 
ing  n(r)  for  all  r  within  rmin  <  r  <  rmMX  at  t  =  0.  The  boundary  condition 
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j  =  0  (26) 

with  j  defined  by  Equation  (4),  corresponds  to  the  case  of  a  perfectly  reflecting  "wall"  in  r  space, 
while  the  condition 


n  =  0  (27) 

corresponds  to  a  perfectly  absorbing  wall.  We  assume  that  no  pores  can  exist  for  r  <  rmjn,  where  rmin 
is  a  parameter  of  the  model.  Packing  constraints  require  rmin  to  be  somewhat  larger  than  the  size  of 
the  hydrophilic  headgroups  (=  0.7  nm)  that  make  up  the  surface  of  the  pore,  since  the  pore  wall  must 
consist  of  at  least  several  phospholipid  molecules25.  For  our  calculations  we  use  rmm  =  1 .0  nm. 

The  problem  of  pore  creation  and  destruction  is  unsolved.  Theories  based  on  treatment  of  insta¬ 
bilities  in  planar  membranes  can  identify  conditions  for  instability  onset26'28,  but  do  not  show  that 
such  an  instability  actually  results  in  pore  structures.  The  formation  of  depressions  or  dimples29  and 
membrane  breathing  modes30,  which  might  precede  pore  structure  formation,  has  been  considered. 

None  of  this  work  is  sufficiently  advanced  to  permit  one  to  compute  pore  birth  or  death  rates.  Wc  are 
therefore  forced  to  make  ad  hoc  assumptions  regarding  the  birth  and  death  rate.  We  assume  that  the 
number  of  pores  changes  due  to  the  creation  and  destruction  of  pores  of  radius  rmm,  with  the  conse¬ 
quence  that  the  boundary  condition  at  r  -  rmin  is 

j  -  K  -N4  at  r  =  rmin  (28) 

where  Ne  and  are  the  creation  and  destruction  rates  of  pores  of  radius  rmjn. 

To  form  a  pore,  which  can  be  thought  of  as  an  excitation  of  the  membrane,  the  membrane  must 
go  through  configurations  which  are  energetically  unfavorable.  We  characterize  this  transition  by  intro¬ 
ducing  a  potential  barrier  A.  Wc  expect  that  this  potential  barrier  will  be  a  function  of  U .  Because 
the  electrical  energy  is  proportional  to  U 2  (for  r  =  re  V  =  V p)  we  further  assume  that 
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A  =  5f  -  a  U1 

where  Sc  and  a  are  constants.  Introducing  an  Arhcnius  factor,  wc  write 


(29) 


’ 

v  exp  - 


(30) 


where  the  attempt  rate  v  is  a  constant  with  dimensions  of  inverse  time.  This  expression  is  valid  as 
long  as  the  fraction  of  the  membrane  area  occupied  by  pores  is  small.  Assuming  that  the  probability 
that  a  pore  of  radius  rmjn  is  destroyed  is  independent  of  Up.  wc  write 


X  n  (rmin)  exp 


(31) 


where  x  is  a  constant  with  dimensions  of  velocity.  Using  Equations  (4),  (29),  (30)  and  (31),  the  boun¬ 
dary  condition  at  r  =  rmin  becomes 


5/  ,  /  _3a£ 

*  * 

ii 

5C  -  aU2 

dr  kT  dr 

=  Z/(rmin)cxp 

kT 

W  J 

-  v  exp 

{  kT  J 

We  now  consider  the  boundary'  condition  at  r  =  rmlx  and  the  related  issue  of  large  pores.  In  the 
follow  paragraphs  we  argue  that  the  large  pores,  whose  unlimited  growth  causes  membrane  rupture, 
form  a  distinct  sub-population  that  must  be  treated  separately  form  the  small  pores,  and  we  derive  an 
equation  to  describe  them.  The  limiting  size  separating  the  large  pores  from  the  small  pores  is  rm„, 
and  the  boundary  condition  that  we  apply  at  rmn  joins  together  the  two  sub-populations. 

3A£ 

We  begin  by  defining  rc (U)  to  be  the  largest  value  of  r  for  which  — - —  =  0.  For  r  >  rc,  AE 

or 

is  a  monotonicallv  decreasing  function  of  r.  This  means  that  a  large  pore  is  unstable,  and  can  grow 
until  the  membrane  ruptures.  For  r  sufficiently  large,  the  "drift  flux"  term  in  Equation  (4)  is  much 
larger  than  the  "diffusion  flux"  term,  and  the  probability  that  such  a  large  pore  will  decrease  in  size  is 
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very  small. 

For  the  description  of  the  small  pores,  it  is  therefore  appropriate  to  use  any  value  of  rmM  that  is 
sufficiently  large,  and  to  apply  the  "absorbing  waH”  condition  (Equation  (27))  there.  In  the  numerical 
simulations,  we  set  rmM  =  2  rc( 0).  The  presence  of  a  pore  of  radius  2rc(0)  (  =  2X1CT6  cm)  docs  not 
mean  that  the  membrane  has  already  ruptured,  as  such  a  pore  is  still  very  much  smaller  than  the  mem¬ 
brane  (=  .1  cm);  rather,  it  implies  that  the  membrane  is  very  likely  to  rupture  soon,  since  the  probabil¬ 
ity  of  such  a  pore  shrinking  again  is  very  small.  Until  the  membrane  docs  rupture,  the  contribution 
from  the  large  pore  or  pores  will  dominate  the  membrane  conductance  G ,  and  therefore  determine  the 
time  evolution  of  U . 

Because  our  model  is  statistical  in  nature  and  a  single  large  pore  can  dominate  G ,  we  must  con¬ 
sider  the  effect  of  fluctuations  in  the  number  of  large  pores.  We  have  defined  n(r)  dr  to  be  the 
number  of  pores  with  radius  between  r  and  r  +  dr.  To  be  more  exact,  n{r)  dr  is  the  average  over  an 
ensemble  of  a  large  number  of  identical  membranes  of  the  number  of  pores  with  radius  between  r  and 
r  +  dr ;  the  actual  number  in  any  particular  membrane  will  in  general  differ  from  n(r)  dr .  If  the 
number  of  pores  is  large  and  the  conductance  of  each  pore  is  small,  the  pore  populations  of  different 
members  of  the  ensemble  are  similar  and  all  members  arc  well  described  by  the  ensemble  average.  On 
the  other  hand,  if  the  number  of  pores  is  small  and  the  conductance  of  each  pore  is  large,  the  pore 
populations  of  different  members  of  the  ensemble  are  quite  different,  and  no  member  of  the  ensemble 
is  well  described  by  the  ensemble  average. 

When  one  or  more  large  pores  dominate  G ,  the  growth  of  the  large  pores  determines  the  time 
dependence  of  G .  For  this  reason  it  is  important  to  consider  the  time  evolution  of  the  radius  of  a  large 
pore.  We  start  by  observing  that  Smoluchowski’s  equation  (Equation  (3))  is  a  limiting  case  of  the 
solution  of  the  problem  of  a  biased  random  walk31.  Consider  a  particle  of  mass  m  that  obeys  the  gen¬ 
eralized  Langevin  equation 
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d2r' 


kT  dr' 
Dp  dt 


dAE 
dr ' 


+  g(t) 


(33) 


where  g(t)  is  a  random  function  that  varies  rapidly  in  time.  Proper  statistical  treatment  of  Equation 
(33)  leads  to  the  Fokker-Planck  equation.  Our  system  corresponds  to  a  "particle”  with  no  inertia,  so 
that  statistical  treatment  leads  to  the  Smoluchowski  equation  (3).  In  the  limit  where  m  lends  to  zero 
and  \dAE!dr’\  »•  |g(/)j.  Equation  (33)  becomes 


dr'  _  DP  dAE  n 

dt  ~  kT  dr'  K 

Thus  we  see  that  "velocity"  of  the  "particle"  is  determined  by  the  balance  between  the  drag  and  the 
driving  force. 

In  the  interest  of  clarity,  wc  also  give  an  alternative  derivation  of  Equation  (34).  In  general,  r\ 
the  expectation  value  of  the  radius  of  a  pore  in  a  membrane  described  by  the  density  function  n ,  is 
defined  by 


j"  r  n(r,t )  dr 

'  K 


where  Np  is  the  number  of  pores  in  the  membrane,  i.e. 


(35) 


Np  ~  J  n(r-‘ 


t )  dr 


(36) 


We  seek  an  equation  for  the  evolution  of  r'  with  time.  Consider  a  membrane  with  a  single  pore  of 
radius  r ,  (>  rc  (0))  at  /  =0.  The  initial  value  of  the  density  function  is 


n  =  6(r  —  r  j)  (37) 

where  8(rj  is  the  Dirac  delta  function.  Note  that  A‘p  =  1  at  t  =  0.  We  assume  that,  throughout  its 
evolution,  both  n  (t)  and  are  zero  at  the  limits  of  the  range  of  r.  To  derive  the  desired  equation. 
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we  take  moments  of  Equation  (3).  We  start  by  integrating  Equation  (3)  over  the  entire  range  of  r. 

dN 

The  left  member  of  the  resulting  equation  is  — — ,  the  time  derivative  of  the  number  of  pores.  When 

dt 

we  integrate  by  parts,  we  see  that  the  right  member  is  zero.  The  resulting  equation  is  easily 
integrated,  yielding 


Np  =  1  (38) 

Next  we  multiply  Equation  (3)  by  r  and  integrate  over  r .  The  left  member  of  Equation  (3)  is  now  just 

the  time  derivative  of  the  expectation  value  of  r.  Integrating  the  right  member  by  parts,  the  first 
dt 

term  vanishes  and  the  second  term  becomes 


f,  A 

dAE 

J  dr 

n  ~ 

dr 


n 

kT  J 


3A E 
dr 


dr 


(39) 


Since  it  is  the  diffusion  term  that  causes  n(r,t )  to  spread,  and  the  diffusion  term  is  small  compared  to 
the  drift  term  for  r  »  /?c <oj.  we  expect  n  to  remain  sharply  peaked  during  its  subsequent  evolution.  In 
dAE 

this  case,  the  term  — —  can  be  removed  from  the  integral  with  the  understanding  that  it  is  to  be 

evaluated  at  r  =  r'.  The  remaining  integral  is  equal  to  unity,  and  the  resulting  equation  is  clearly 
equivalent  to  Equation  (34). 

We  now  solve  Equation  (34).  Substitution  of  Equation  (6)  for  A E  into  Equation  (34)  and  use  of 
Equations  (7)  -  (24)  yields  an  equation  that  depends  upon  U,  and  therefore  is  coupled  to  the  equation 
that  describes  the  electrical  circuit.  But  for  sufficiently  large  pores  A Ee,  which  contains  the  U- 
dependcnce,  is  small  compared  to  A £w  (compare  Equation  (7)  to  Equations  (22)  and  (23)),  and  can  be 
neglected.  The  resulting  equation  is  independent  of  U  and  can  be  solved  analytically.  Setting  A E 
equal  to  A EM  and  using  Equation  (7)  transforms  Equation  (34)  into 


dr' 

dt 


(Y-ro 
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Equation  (40)  is  separable,  and  its  solution  is 


r'  =  r£(0)  +  (ri  -  rc(0»  exp(— ) 

xd 


where  the  characteristic  time  is 


(41) 


kT 

2  nTDp 


(42) 


and  rt  is  the  initial  value  of  r\ 

Note  that  Equation  (41),  which  describes  the  evolution  of  large  pores,  is  deterministic.  To  com¬ 
bine  the  statistical  description  of  the  small  pores  with  the  deterministic  description  of  large  pores,  we 
need  a  criterion  for  the  time  of  formation  of  large  pores.  The  probability  that  a  pore  becomes  large 
during  the  interval  between  t  and  ;  +  dt  is  j(rmtx,t)dt,  where  j  is  defined  by  Equation  (4).  Since 
n  =  0  at  r  =  only  the  first  term  in  Equation  (4)  contributes.  Since  large  pores  cannot  be  des¬ 
troyed  (unless  the  membrane  ruptures)  and  the  probability  of  an  additional  large  pore  forming  does  not 
depend  upon  how  many  large  pores  already  exist,  the  formation  of  large  pores  is  a  Poisson  process. 
The  probability  P(Nh  t)  that  N,  large  pores  exist  at  time  t  is  given  by  the  Poisson  distribution 


where  the  parameter  p  is 


P{Nt,t)  =  e* 


(43) 


4  = 


dt' 


The  probability  that  the  membrane  has  no  large  pores  at  lime  t  is  simply 


P(0,t)  =  exp(-p) 


(44) 


(45) 


Anril  5.  19SQ 


Page  18 


The  dominant  influence  of  a  small  numbers  of  large  pores  on  the  rupture  process  has  an  impor¬ 
tant  consequence;  thermodynamic  fluctuations  cause  the  lifetime  of  a  membrane  undergoing  rupture  to 
vary  from  experiment  to  experiment,  even  if  the  experimental  setups  are  identical.  This  is  due  to  the 
fact  that  the  creation  of  the  large  pore  that  causes  rupture  is  a  random  process.  The  straightforward 
way  to  compute  r|  t  the  creation  time  of  the  i-th  large  pore,  is  to  use  a  random  number  generator.  For 
the  purposes  of  this  paper,  it  is  more  convenient  to  use  a  completely  deterministic  model.  We  therefore 
require  r,  to  equal  that  value  of  t  for  which  p  equals  »;  that  is,  the  first  large  pore  is  created  at  rj, 
when  p  =  1,  the  second  at  f2.  when  p  =  2,  etc,  where  p  is  defined  by  Equation  (44).  This  is  a  reason¬ 
able  scheme,  since  the  Poisson  distribution  P(Nlt  t)  (Equation  (43))  has  the  property  that  p  is  the 
expectation  value  of  A', ,  the  number  of  large  pores. 

Once  a  large  pore  forms,  its  radius  as  a  function  of  time  is  given  by  Equation  (41).  We  therefore 
define  the  density  function  n‘  as 


n  -  n  -r  £  5(r  -  r\i  -  /,))  (46) 

l 

where  n  obeys  Equation  (3)  for  rmjn  <  r  <  rmv,  and  equals  zero  for  r  >  rm„.  and  r,  is  the  creation 
time  of  the  i  -th  large  pore. 

This  division  of  n  into  two  pieces  has  one  drawback;  if  the  conductance  of  a  single  pore  of 
radius  rmtx  is  a  sizable  fraction  of  the  conductance  of  the  membrane,  there  will  be  a  noticeable  discon¬ 
tinuity  in  C  (?)  at  t  =  tj.  The  only  way  to  avoid  this  discontinuity  and  still  correctly  model  the  effect  of 
the  unstable  large  pores  (a  pore  with  radius  sufficiently  large  will  always  dominate  G )  is  to  reformu¬ 
late  the  problem  and  perform  a  Monte  Carlo  calculation,  which  requires  much  more  computing  power. 
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For  the  initial  condition,  we  assume  that  at  t  =  0  the  membrane  is  in  equilibrium  with  U  =  0. 
The  quasi-steady  state  solution  to  Equation  (3)  can  be  found  analytically;  it  is 


n(r)  = 


—exp 

X 


hd  -5C  -  2 Jttfr  -  rmin)  +  nT(r2  -  rmin2) 


kT 


J  ''min  £  r  <  rc(0) 
r  >  rc  (0) 


(47) 


THE  EXTERNAL  CIRCUIT 

To  close  the  system  of  equations,  one  needs  an  additional  equation  to  describe  the  time  evolution 
of  U .  The  membrane  is  adequately  modeled  as  a  capacitor  C  in  parallel  with  a  resistor  with  conduc¬ 
tance  G(t).  The  value  of  C  can  be  computed  with  the  familiar  formula  from  elementary  electrostatics 


C 


4  r*,Am 
h 


(48) 


where  Am  is  the  area  of  the  membrane.  The  capacitance  of  the  partition  that  holds  the  membrane  is 
much  smaller  than  the  capacitance  of  the  membrane  itself  and  can  be  neglected.  G(r)  can  be  expressed 
in  terms  of  an  integral  over  the  pore  population.  The  current  lp  flowing  through  a  pore  of  radius  r  in 
a  membrane  with  voltage  U  across  it  is 


/  _ _ SL _ 

P  RAr)  +  Rp(r) 


(49) 


where  Rs(r )  and  Rp(r)  arc  defined  in  Equations  (19)  and  (20),  respectively.  The  current  /  flowing 
through  the  entire  membrane  is  simply  the  sum  of  the  current  flowing  through  each  pore  and  can  be 
expressed  as  an  integral  over  n(r,  r) 


/  =  U 


s 


n  (r.t) 


t-l 


Rf(r)  +  Rp(r) 


dr 


The  membrane  conductance  G  -  -jj  is  therefore 


(50; 
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C(t) 


■J 


n  (r  ,t) 


Rs(r)  +  Rp(r) 


-1 


dr 


(51) 


Figure  3  is  a  schematic  diagram  of  the  entire  circuit.  A  square  wave  pulse  of  amplitude  /,  is 
applied  to  the  electrodes  for  0  <  t  <  tpulse ,  and  at  t  =  tpulse  the  switch  is  opened.  Thus,  charge  is 
supplied  to  the  membrane  from  an  external  source,  but  discharge  after  the  switch  opens  at  the  end  of 
the  pulse  can  occur  only  through  the  membrane.  Such  a  circuit  is  described  by  the  differential  equa¬ 
tions 


C 


dU 

dt 


1P  rn 
Re  +  Rs 


-  U  (G  (r )  + 


1 

Re  +  Rn 


) 


u_ 

R 


if  f  <  tpulse 

if  I  >  tpulse 


and  the  initial  condition 


U(0)  =  0 


NUMERICAL  SOLUTIONS  OF  THE  EQUATIONS 


(52) 


(53) 


In  the  present  model,  the  membrane  is  characterised  by  eleven  parameters  {h ,  e;.  A,  y,  T,  D ,  v, 

X.  5C.  8*.  and  r^).  Unfortunately,  it  is  difficult  to  obtain  better  than  an  order  of  magnitude  estimate 
for  eight  of  them  from  measurements  of  membrane  properties  other  than  electroporation.  In  principle, 
one  could  determine  their  values  by  obtaining  a  "best  fit"  to  experimental  data.  Because  the  number  of 
data  points  one  can  obtain  (the  membrane  potential  U(t )  for  many  values  of  t  for  different  pulse 
lengths,  applied  potentials,  solutes,  temperatures,  etc.)  greatly  exceeds  eight,  the  system  is  multiply 
over  determined.  Therefore,  unless  the  system  is  mathematically  ill-conditioned,  a  fitting  procedure 
could  be  used  to  determine  the  values  of  the  parameters.  Due  to  the  large  amount  of  computing  time 
required  for  such  a  fitting  proceedurc,  we  do  not  fit  data.  Instead,  we  show  that,  for  one  plausible  set 
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of  values  of  the  parameters  (the  "standard  membrane"),  the  present  model  describes  both  REB  and  rup¬ 
ture.  In  future  work  we  will  explore  the  behavior  of  the  solutions  as  the  parameters  are  varied. 

We  solve  the  Equation  (3)  using  the  Crank-Nicholson  method.  In  addition  to  the  eleven  parame¬ 
ters  that  describe  the  membrane,  one  must  describe  the  remainder  of  the  experimental  apparatus  by 
specifying  the  temperature  T,  the  electrical  conductivity  of  the  solution  o,  the  radii  r;  and  charge  of 
each  species  of  small  ion,  the  shunt  resistance  RN ,  the  resistance  RE  of  the  electrodes  and  the  solution, 
and  the  strength  of  the  current  source  as  a  function  of  time.  We  use  a  square  wave  of  amplitude  /, 
and  duration  tpulse  for  the  applied  current.  We  specify  the  grid  size  A r  =  (rmiX  -  rmin)/m  and  the 
time  step  size  At  =  Ar2  dt!Dp  by  quoting  the  values  of  the  parameters  m  (the  number  of  intervals 
between  grid  points)  and  dt  (the  time  step  size  in  units  of  the  diffusion  time  based  on  the  grid  spac¬ 
ing).  The  values  of  all  of  the  parameters  used  in  the  numerical  calculations  are  given  in  Table  1. 

RESULTS 

Figures  5-12  show  the  results  of  our  numerical  calculations.  We  simulated  five  charge  pulse 
experiments,  which  differed  only  in  the  amount  of  injected  charge  Q  =  tpuiseIo-  The  pulse  length  ipuist 
was  0.4jis  for  all  five  simulations,  the  amount  of  injected  charge  Q  was  5  nanocoulombs,  10  nano¬ 
coulombs,  15  nanocoulombs,  20  nanocoulombs,  and  25  nanocoulombs,  and  the  membrane  parameters 
are  given  in  Table  1.  Figures  5  and  6  snow  U (r ),  Figures  7  and  8  show  ,V(r),  and  Figures  9  and  10 
show  G(:)  for  all  five  simulations.  The  curves  in  Figures  5  -  10  are  labeled  by  values  of  the  injected 
charge.  There  are  two  time  scales  of  interest:  a  short  lime  scale  fapprox  1  us  )  characteristic  of  mem¬ 
brane  charging  and  REB,  and  a  long  time  scale  (=80  jjls)  characteristic  of  membrane  recovery  or  rup¬ 
ture.  Figures  5.  7  and  9  show  short  time  scale  behavior,  while  Figures  6,  8,  and  10  show  long  time 
scaic  behavior. 

We  denote  by  l\,  the  value  of  the  membrane  potential  at  the  end  of  the  applied  pulse  i:  =  0.4p.o. 
Viewed  as  a  function  of  Q,  i'0  increases  for  small  values  of  Q ,  reaches  a  maximium  value,  and  then 
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decreases  again.  We  denote  by  U0,c  the  maximum  value  of  Uo,  and  by  Qe  the  maximum  value  of  Q. 
For  these  simulation,  Qc  =  20  nanocoulombs  and  i/0,e  “  0-94  v.  If  the  injected  charge  is  greater  than 
Qc ,  REB  occurs  before  the  end  of  the  pulse,  and  the  membrane  conductivity  increases  so  much  that 
the  membrane  starts  to  discharge  before  the  pulse  ends  (see  Figs.  5,  7,  and  9).  A  membrane  undergo¬ 
ing  REB  discharges  completely  in  less  than  1  ps.  Following  discharge  the  membrane  recovers  in 
about  80  ps.  The  recovery'  can  best  be  seen  in  Figures  8  and  10.  The  number  of  pores  increases  in 
less  than  1  ps,  which  looks  like  a  step  function  in  the  long  time  scale  plot  (Fig.  8),  and  then  decays 
exponentially  back  to  its  equilibrium  value.  The  membrane  conductance  also  increases  in  about  1  ps, 
but  it  decays  in  two  stages,  as  can  be  seen  in  Figure  10.  The  first  stage,  during  which  the  decay  is 
rapid,  is  caused  by  the  shrinking  of  the  the  pores;  the  second  stage,  during  which  the  decay  is  much 
slower,  is  caused  by  the  decrease  in  the  number  of  pores.  The  membrane  returns  to  its  initial  state  in 
about  80  ps.  The  evolution  of  G  and  A’  during  REB  are  best  understood  by  considering  the  evolution 
of  n(r)  shown  in  Figure  11.  Figure  1 1  shows  n(r)  at  different  times  between  0  and  60  ps  for  Q  -  20 
nanocoulombs;  the  curves  in  Figure  1 1  are  labeled  with  the  time  since  the  begining  of  the  pulse. 

Driven  by  the  strong  electrical  forces,  n  increases  rapidly  for  all  r  until  REB  occurs.  As  the  mem- 

brane  begins  to  discharge,  the  pore  creation  rate  drops,  but  the  force - - —  still  acts  to  increase  the 

or 

size  of  existing  pores  until  r  =  0.5  -  0.6  ps,  by  which  time  the  membrane  has  discharged  sufficiently 
that  the  electrical  force  no  longer  dominates  the  mechanical  force,  which  favors  contraction.  The  pores 
shrink  faster  than  they  are  destroyed,  causing  a  decrease  in  the  average  pore  size.  During  this  lime,  N 
decreases  slowly  while  C  decreases  rapidly  due  to  the  rapid  shrinking  of  the  pores.  Eventually  n 
relaxes  to  a  quasi-static  equilibrium;  its  shape  no  longer  changes,  but  its  magnitude  decreases  due  to 
destruction  of  pores  radius  rmin.  During  this  final  phase.  G  decrease  slowly,  its  variation  entirely  due 
to  the  change  in  A’,  as  the  membrane  asymptotically  approaches  its  original  state. 

For  injected  charge  less  than  Qc ,  REB  docs  not  occur  during  the  pulse,  and  the  membrane  retains 
its  charge  for  a  much  longer  time.  In  this  case,  the  number  of  pores  and  the  membrane  conductance 
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continue  to  grow  even  after  the  end  of  the  pulse.  There  even  exists  a  range  of  Q  for  which  the  mem¬ 
brane  does  not  undergo  REB  until  after  the  end  of  the  pulse.  For  Q  =  15  nanocoulombs,  the  mem¬ 
brane  discharges  rapidly  for  t  between  0.8  ps  and  4  ps,  at  which  time  U  =  0.05  v  and  the  membrane 
conductance  only  slightly  exceeds  the  conductance  of  an  unexcited  membrane  (=  4.5  x  10"9).  The 
membrane  eventually  discharges  completely  with  an  RC  time  constant  of  about  2  s. 

For  Q  sufficiently  small,  no  breakdown  occurs.  For  (2=5  nanocoulombs,  the  membrane  rapidly 
charges  up  to  U  =  0.25  v.  At  f  =  80  ps,  the  membrane  pore  population  is  approaching  a  quasi-steady 
state,  and  the  conductance  has  increased  by  a  factor  of  about  10.  The  membrane  then  starts  to 
discharge  with  a  time  constant  of  about  0.1  s.  As  the  membrane  discharges,  G  decreases  and  the  rate 
of  discharge  slows  until  it  approaches  the  unexcited  membrane  RC  time  constant  of  about  2  s. 

Intermediate  values  of  Q  lead  to  the  far  more  interesting  case  of  mechanical  rupture.  For  Q  = 

10  nanocoulombs,  the  membrane  charges  up  to  U  =  0.5  v.  It  remains  almost  constant  for  about  5  ps, 
during  which  time  the  membrane  conductance  has  increased  sufficiently  to  cause  the  membrane  to 
discharge.  As  the  membrane  discharges,  G  reaches  a  maximum  and  begins  to  decrease  at  about  r  = 

15  ps.  Until  t  =  30  ps  the  evolution  resembles  the  case  Q  =  15  nanocoulombs,  slowed  down  by  a 
factor  of  about  5.  The  resemblance  ceases  at  t  =31  ps  at  which  time  G  starts  to  increase  again  due  to 
the  effect  of  an  unstable  pore.  (The  discontinuities  in  the  curve  in  Figure  10  are  artifacts  of  the  tech¬ 
nique  used  to  keep  track  of  the  large  pores;  in  a  real  membrane  G(r)  is  continuous.)  As  G  increases 
the  membrane  discharges  more  and  more  rapidly,  until  it  is  completely  discharged  by  t  =  72  ps.  The 
unstable  pore  continues  to  grow  until  the  membrane  ruptures.  The  time  constant  for  the  exponential 
growth  of  the  unstable  pore  is  given  by  Equation  (42).  For  the  present  example,  %  -  13  ps.  The 
evolution  of  n  leading  to  mechanical  rupture  is  shown  in  Figure  12. 

DISCUSSION 
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We  have  presented  a  quantitative  theory  of  electroporation  of  lipid  bilayer  membranes.  The 
theory  successfully  describes  the  four  different  fates  of  membranes  in  charge  pulse  experiments.  We 
have  not  fit  any  data,  but  rather  we  have  performed  numerical  "experiments"  with  a  "standard  mem¬ 
brane".  While  our  results  qualitatively  agree  with  the  experimental  work  of  Benz  et  al\  comparison  of 
our  results  with  theirs  shows  the  following  differences: 

(1)  The  membrane  voltage  at  the  end  of  the  pulse  U  ( tpulu)  for  our  "standard  membrane"  is  a  much 
stronger  function  of  the  injected  charge  than  is  observed  experimentally. 

(2)  The  time  delay  before  membrane  rupture  for  our  "standard  membrane"  is  shoner  by  a  factor  of 
about  5  than  is  observed  experimentally. 

We  do  not  know  how  much  the  agreement  with  experiment  can  be  improved  by  changing  the  values  of 
the  parameters;  we  plan  to  explore  this  question  in  the  future. 

If  the  data  cannot  be  fit  better  by  changing  the  values  of  the  parameters,  there  arc  several  areas  in 
which  the  model  can  be  improved.  Improved  approximations  for  AE  and  better  modeling  of  pore  crea¬ 
tion  and  destruction  might  improve  agreement  with  experiment.  In  addition,  the  current-voltage  rela¬ 
tion  of  a  pore  is  probably  not  linear12’32,  and  this  might  effect  the  results. 

Another  topic  for  future  work  is  modeling  voltage  clamp  experiments,  which  deal  with  smaller 
trans-membrane  voltages  and  longer  time  scales.  Unfortunately,  simulation  of  a  voltage  clamp  experi¬ 
ment  with  our  present  computer  program  takes  a  prohibitively  long  time;  we  would  need  to  reduce  the 
running  time  by  a  factor  of  1000. 

In  summary,  we  believe  that  the  present  theory  correctly  describes  the  overall  features  of  electro¬ 
poration.  but  further  work  is  required  to  obtain  detailed  quantitative  agreement  with  experiment. 
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|  Table  1 

Am 

membrane  area 

C 

capacitance  of  membrane 

9.61  nF 

D 

diffusion  constant 

dt 

time  step  size  (in  units  of  £>/(Ar)2) 

0.5 

h 

membrane  thickness 

2.8  nm 

kT 

temperature 

IlllliiCfSi 

Em 

series  resistance  of  electrolyte,  electrodes,  and  wires 

30  n 

mm 

internal  resistance  of  current  source 

50  a 

El 

large  pore  creation  size 

40  nm 

minimum  pore  radius 

1  nm 

r+ 

radius  of  positive  ions 

0.2  nm 

r_ 

radius  of  negative  ions 

0.2  nm 

t pulse 

source  current  pulse  length 

0.4  ps 

z* 

charge  of  positive  ions  (in  units  of  proton  charge) 

+1 

z_ 

charge  of  negative  ions  (in  units  of  proton  charge) 

-1 

nH 

pore  edge  energy  density 

2  pergs/cm 

r 

membrane  surface  tension 

5e 

pore  creation  energy  barrier 

2.04  10-12  ergs 

!  §4 

pore  destruction  energy  barrier 

2.04  10~12  ergs 

A  r 

grid  spacing 

.0195  nm 

£/ 

dielectric  constant  of  lipid 

2.1  Eo 

dielectric  constant  of  water 

80  Eo 

v 

pore  creation  rate  prefactor 

1028  s"1 

a. 

conductivity  of  bulk  solution 

.098  Q_1  cm-1 

Lz _ 

pore  destruction  rate  prefactor 

1018  cm/s 
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List  of  Symbols  (see  also  Table  1) 
coefficient  of  U 2  in  A 
concentration  of  i-th  ion  species 
proton  charge 
electric  field 

electric  field  inside  a  pore 
membrane  conductance 

steric  hindrance  factor  for  ion  of  radius  r,  in  pore  of  radius  r 

applied  current 

measured  current 

current  through  a  pore 

flux  density  of  pores  in  radius  space 

probability  density  function  of  pores  (small  pores  only) 

probability  density  function  of  pores  (including  large  pores) 

number  of  pores  in  the  membrane  as  a  function  of  time 

pore  creation  rate 

pore  destruction  rale 

Parsegian’s  function  for  Bom  energy 

pore  radius 
radius  of  ion  of  type  i 
pore  resistance 

spreading  resistance  (includes  both  sides  of  the  membrane) 

time  since  start  of  pulse 

transmembrane  potential 

applied  voltage 

measured  voltage 

transmembrane  potential  at  the  end  of  a  square  pulse 
maximum  transmembrane  potential  at  the  end  of  a  square  pulse 
voltage  drop  across  a  conducting  pore 

voltage  drop  in  electrolyte  near  the  ends  of  a  conducting  pore  (both  ends) 
pore  "energy" 

difference  in  electrical  pressure  at  the  pore  edge 
electrical  mobility  of  the  i-th  ionic  species 
potential  barrier  for  pore  formation 

standard  chemical  potential  of  an  ion  of  species  i  inside  a  pore 

electrical  conductivity 

conductivity  of  electrolyte  within  a  pore 

time  constant  for  growth  of  an  unsable  pore 
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FIGURE  CAPTIONS 

FIGURE  1.  The  experimental  apparatus  consists  of  a  vessel  filled  with  an  electrolytic  solution. 

The  vessel  is  separated  into  two  compartments  by  a  partition.  On  each  side  of  the  partition  is  a  planar 
electrode  connected  to  a  power  supply  and  a  measuring  circuit.  The  membrane  spans  an  aperture  in 
the  partition. 

FIGURE  2.  An  equivalent  circuit  for  a  voltage  clamp  experiment  using  the  apparatus  shown  in 
Figure  1.  The  voltage  source  f/,(r)  in  series  with  the  resistor  RT  models  the  power  supply,  the  resistor 

Re  models  the  electrolyte,  electrodes  and  wires,  and  the  capacitor  C  in  parallel  with  the  resistor  ^ 

model  the  membrane.  The  ammeter  measures  the  current  /m(r)  flowing  through  the  electrodes. 

FIGURE  3.  An  equivalent  circuit  for  a  charge  pulse  experiment  using  the  apparatus  shown  in  Fig¬ 
ure  1.  The  current  source  /,(r)  in  parallel  with  the  resistor  R#  models  the  power  supply,  the  resistor 

R£  models  the  electrolvie,  electrodes  and  wires,  and  the  capacitor  C  in  parallel  with  the  resistor  — 

G(:) 

model  the  membrane.  The  switch  is  opened  at  the  end  of  the  current  pulse.  The  voltmeter  measures 
the  potential  difference  Um(t)  between  the  electrodes. 

FIGURE  4.  Pore  "energy"  vs.  radius  for  a  membrane  described  by  the  parameters  listed  in  Table  1. 
Each  curve  is  labeled  by  the  corresponding  membrane  voltage,  which  ranges  between  0  and  1  volt.  A 
pore  of  radius  rmin  is  stable  only  for  U  less  than  about  0.35  v,  and  the  force  favors  expansion  of  pores 
larger  than  20  nm  for  all  values  of  U . 

FIGURE  5.  Membrane  voltage  vs  time  (short  time  scale),  for  a  simulated  charge  pulse  experiment. 
Each  curve  is  labeled  by  the  corresponding  value  of  the  injected  charge  Q .  The  curves  for  Q  =  25 
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nanocoulombs  and  20  nanocoulombs  show  REB  while  the  others  do  not  The  values  of  the  membrane 
parameters  are  given  in  Table  1. 

FIGURE  6.  Membrane  voltage  vs  time  (long  time  scale)  for  a  simulated  charge  pulse  experiment. 
Each  curve  is  labeled  by  the  corresponding  value  of  the  injected  charge  Q .  The  curves  for  Q  =25 
nanocoulombs  and  20  nanocoulombs  arc  the  spikes  at  t  =  0.  The  curve  for  Q  =  15  nanocoulombs 
shows  that  the  membrane  suffered  REB  at  r  =  2ps,  but  the  membrane  recovered  before  it  had  time  to 
discharge  completely.  The  curve  for  Q  =  10  nanocoulombs  shows  rupture,  while  the  curve  for  Q  -5 
nanocoulombs  shows  that  the  membrane  conductance  did  not  increase  enough  to  discharge  the  mem¬ 
brane.  The  values  of  the  membrane  parameters  are  given  in  Table  1. 

FIGURE  7.  Number  of  pores  vs  time  (short  time  scale)  for  a  simulated  charge  pulse  experiment 
Each  curve  is  labeled  by  ihe  corresponding  value  of  the  injected  charge  Q .  For  Q  -  25  nanocoulombs, 
,Y  increases  from  about  6  to  about  2  x  108.  The  values  of  the  membrane  parameters  are  given  in  Table 
1. 


FIGURE  8.  Number  of  pores  vs  time  (long  time  scale)  for  a  simulated  charge  pulse  experiment. 
Each  curve  is  labeled  by  the  corresponding  value  of  the  injected  charge  Q .  For  Q  -25  nanocoulombs 
and  20  nanocoulombs,  cases  for  which  REB  occurs,  N  increases  to  about  10s  in  less  than  .5  ps  and 
then  decays  exponentially  with  a  time  constant  of  =  4.5  (is.  For  Q  -  15  nanocoulombs,  N  increases 
rapidly  to  about  10s  and  remains  almost  constant  for  about  4  ps  before  the  exponential  decrease.  For 
Af  =  10  nanocoulombs.  N  increases  to  about  2  x  103  in  about  5  p  and  remains  almost  constant  for 
about  30  ps  before  the  decay  phase.  The  membrane  in  this  case  ruptures.  For  (2=5  nanocoulombs, 
S  increases  to  about  40  in  80  ps.  .V  will  return  to  its  initial  value  as  the  membrane  discharges  with  a 
time  constant  of  about  2  s.  The  values  of  the  membrane  parameters  are  given  in  Table  1. 
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FIGURE  9.  Membrane  conductance  vs  time  (short  time  scale)  for  a  simulated  charge  pulse  exper¬ 
iment.  Each  curve  is  labeled  by  the  corresponding  value  of  the  injected  charge  Q .  For  membranes  that 
undergo  REB  ( Q  -  25  nanocoulombs,  20  nanocoulombs),  G  increases  by  8-9  orders  of  magnitude  in 
less  that  .3  ps.  The  values  of  the  membrane  parameters  arc  given  in  Tabic  1. 

FIGURE  10.  Membrane  conductance  vs  time  (long  time  scale)  for  a  simulated  charge  pulse  experi¬ 
ment.  Each  trace  is  labeled  by  the  corresponding  value  of  the  injected  charge  Q .  The  curves  for  Q  = 
25  nanocoulombs  and  Q  -  20  nanocoulombs,  the  two  cases  of  REB,  show  a  rapid  rise  in  conductance, 
followed  by  relaxation  to  the  initial  state  in  two  stages.  The  first  stage,  during  which  G  decays 
rapidly,  lasts  about  1.5  ps  and  is  due  to  the  shrinking  of  the  pores.  The  second  stage,  during  which  G 
decays  much  more  slowly  and  approaches  its  initial  value,  lasts  about  80  ps  and  is  due  to  the  death  of 
the  pores.  The  curve  for  Q  -  15  nanocoulombs  is  similar  to  the  higher  voltage  curves,  except  that 
REB  is  delayed  and  occurs  after  the  end  of  the  pulse,  causing  the  first  stage  to  take  7  ps.  The  curve 
for  Q  =10  nanocoulombs  shows  membrane  rupture;  the  conductance  increases  for  i  >  30  ps  due  to 
the  instability  of  a  large  pore.  (The  discontinuities  in  the  curve  are  an  artifact  of  the  algorithm  used  to 
treat  the  large  pores.)  For  Q  -  5  nanocoulombs,  the  conductance  increases  to  a  new  equilibrium  value 
appropriate  to  its  non-zero  membrane  potential.  The  membrane  will  discharge  with  an  RC  time  con¬ 
stant  of  about  2  s  as  the  membrane  conductance  returns  to  its  initial  value.  The  values  of  the  mem¬ 
brane  parameters  are  given  in  Table  1 . 

FIGURE  II.  n  vs.  r  for  a  pore  undergoing  REB.  The  density  function  n  is  of  fundamental  impor¬ 
tance  in  the  present  theory'  of  electroporation,  even  though  it  impossible  to  measure  it  directly.  The 
curves  arc  labeled  by  the  time  in  ps  since  the  start  of  the  charging  pulse  and  represent  n(r.  t)  for  Q  - 
20  nanocoulombs.  The  population  of  pores  of  all  sizes  grows  for  the  duration  of  the  pulse  (which 
corresponds  to  the  time  of  REB).  As  the  membrane  rapidly  discharges  following  REB,  the  pore  des- 
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truction  rate  increase,  causing  a  depletion  in  the  population  of  small  pores.  Large  pores  then  shrink  as 
the  population  relaxes  to  its  initial  state.  The  values  of  the  membrane  parameters  are  given  in  Table  1. 

FIGURE  12.  n  vs.  r  for  a  pore  that  ruptures.  The  curves  arc  labeled  by  the  time  in  ps  since  the  start 
of  the  charging  pulse  and  represent  n(r ,  t)  for  Q  =  10  nanocoulombs.  Note  the  presence  of  a  flux  of 
pores  ( n  has  a  negative  slope)  through  the  boundary  at  r  =  40  nm  for  t  >  30  ps.  This  flux  leads  to 
the  creation  of  the  large  unstable  pores  that  cause  rupture.  The  probability  that  a  large  pore  had  formed 
was  0.63  at  t  =  32  ps.  The  values  of  the  membrane  parameters  arc  given  in  Table  1. 
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